Milestone 1¶
Overview¶
Important Questions¶
Context
Why is this problem important to solve?
Climate change's threat looms ever nearer the entire globe. According to NASA, climate change will cause flooding, drought, heat waves, stronger hurricanes, and a higher sea level in the US. While some climate changes occur naturally, humanity's CO2 emissions already caused climate change that would not have happend otherwise. Susan Solomon et al states the concentration of CO2 in our atmosphere is already so large that the effects will be irreversable for 1,000 years even if all emissions were to stop today. In order to continue enjoying this beautiful planet in the long run, humanity must manage its resources more efficiently.
Objective
What is the intended goal?
The US Environmental Protection Agency estimates that approximately 25% of all greenhouse gas emissions (including CO2) come from electricity production. The same EPA article states that the US produces 62% of its energy by burning fossil fuels such as coal and natural gas. Predicting future CO2 emissions from different energy sources may help policy makers and innovators understand where to focus their green energy efforts.
Key questions
What are the key questions that need to be answered?
- What energy sources emit the most CO2?
- What energy sources emit the least CO2?
- Which energy sources are expected to emit CO2 at an even faster rate in the near future?
- Which energy sources should emit less CO2 as time goes on without intervention?
Problem Formulation
What is it that we are trying to solve using data science?
The purpose of this project is to identify potential opportunities to prevent the acceleration of CO2 emission rates from electricity generation.
Data Dictionary¶
This datset is the past monthly data of Carbon dioxide emissions from electricity generation from the US Energy Information Administration categorized by fuel type such as Coal, Natural gas etc.
- MSN - Reference to Mnemonic Series Names (U.S. Energy Information Administration Nomenclature)
- YYYYMM - The month of the year on which these emissions were observed
- Value - Amount of CO2 Emissions in Million Metric Tons of Carbon Dioxide
- Description - Different category of electricity production through which carbon is emissioned.
Data Preparation¶
Before we visualize and model our data, we must make sure there are no null or otherwise invalid arguments and format the data into a dataframe.
In [18]:
# Set up Google Colab.
from google.colab import drive
drive.mount('/content/drive')
# Avoid scroll-in-the-scroll for large outputs.
from IPython.display import Javascript
def resize_colab_cell():
display(Javascript('google.colab.output.setIframeHeight(0, true, {maxHeight: 5000})'))
get_ipython().events.register('pre_run_cell', resize_colab_cell)
In [19]:
# Import basic libraries.
!pip install statsmodels --upgrade
import pandas as pd
import warnings
import itertools
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from statsmodels.graphics import tsaplots
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from statsmodels.tsa.stattools import coint, adfuller
from statsmodels.tsa.ar_model import AutoReg
warnings.filterwarnings("ignore")
In [20]:
# Load the data.
df = pd.read_excel('/content/drive/MyDrive/Colab Notebooks/MER_T12_06.xlsx')
# Convert dates into standard datetime
df['YYYYMM'] = pd.to_datetime(df['YYYYMM'], errors = 'coerce', format = '%Y%m')
# Drop rows with NaT in date column and set date as index.
df = df.dropna(subset = ['YYYYMM']).set_index('YYYYMM')
df.head(15)
Out[20]:
In [21]:
# Check the info of the dataframe.
df.info()
Observations:
- Our dataframe has 3 object columns with no missing data.
- There are 4,707 observations even after deleting observations with invalid dates.
- Since the Value column is supposed to be a float, we need to convert numeric data in this column to a float and delete all rows with invalid values that cannot be converted to a float.
In [22]:
# Convert valid values in the Value column to float and drop any rows with invalid values in said column.
# Check the info of the dataframe.
df['Value'] = pd.to_numeric(df['Value'], errors = 'coerce', downcast = 'float')
df = df.dropna(subset = ['Value'])
df.head(15)
Out[22]:
Observations:
- Our dataframe now has 2 object columns and 1 float32 column.
- There are now 4,323 observations with no missing data.
- We are now ready to visualize our data!
Dataset Visualization¶
Further Formatting¶
Since the dataset has 8 different sources of carbon emissions, we first need to group by the description. Then we will be able to compare the sources to each other in our visualizatiuons. Once we do so, we can drop the MSN column from the dataset.
In [23]:
# Create a dataframe that groups by description and uses the description as the column names.
EnergyTypes = ['Coal', 'Distillate Fuel', 'Geothermal Energy', 'Natural Gas', 'Non-Biomass Waste',
'Petroleum Coke', 'Petroleum', 'Residual Fuel Oil', 'Total Emissions']
Emissions = pd.DataFrame(index = df.index.unique(), columns = EnergyTypes)
grouped = df.groupby('Description')
for g, c in zip(grouped.groups, EnergyTypes):
Emissions[c] = grouped.get_group(g).iloc[:, 1]
Emissions.head(15)
Out[23]:
Observations:
- The Geothermal Energy and Non-Biomass Waste columns don't have any observations until 1989.
- We will have to drop all rows that have any nulls in order for comparisons to be over the same timeframes and make sense.
In [24]:
# Drop rows with null values in Emissions.
EmissionsDropped = Emissions.dropna()
# View information on our new dataframe.
EmissionsDropped.info()
Observations:
- We now have 331 months of observations for each of the 8 energy types as well as total emissions.
- There are no nulls in our dataset.
In [25]:
# View the head of our data.
EmissionsDropped.head(15)
Out[25]:
Observations:
- Coal often contributes a whole order of magnitude more carbon emissions than all the other energy sources combined.
- Petroleum, Residual Fuel Oil, and Natural Gas all contribute the most after coal to carbon emissions consistently; often an order of magnitude more.
- All of the energy types seem to have peak emissions around July/August and then start to decrease. Later when we check for stationarity, we should check and make sure this isn't a seasonal pattern.
Chart Time¶
Now that the dataframe is in a more user-friendly format, we are ready to create the different plots and graphs.
Box Plots¶
In [26]:
# Create boxplots to visualize each energy source's emissions over this timeframe.
EmissionsDropped.boxplot(figsize=(24, 12), grid = False, fontsize=12, color = 'orange');
Observations:
- Coal has the widest distribution - therefore the greatest standard deviation - of all the energy types.
- Coal is slightly left-skewed.
- Distillate Fuyel, Natural Gas, Petroleum Coke, Petroleum, and Residual Fuel Oil are all slighlty right-skewed.
- Coal is by far the greatest contributor to carbon emissions.
Bar Chart¶
In [27]:
# Compare total emissions from each energy source in a bar plot.
EmissionsDropped.sum().plot.bar(figsize=(24, 12), fontsize=12, ylabel = 'CO2 Emissions (MMT)', xlabel = 'Energy Type', title = '1989 - 2016 CO2 Emissions by Energy Type');
Observations:
- Coal alone contriputes over 3/4 of total CO2 emissions for this time period. Either we heavily rely on coal to generate energy and/or it is highly inefficient.
- Natural gas is the 2nd-greatest contributor. This makes sense given that coal is often used for industrial purposes and most homes and apartments use natural gas for heat and cooking.
- Geothermal energy emits the least amount of CO2. This is probably because geothermal energy is a newer technology, but it may also get more energy per MMT of CO2 emissions.
Line Graph of Each Source Over Time¶
In [28]:
# Show the relationship between power generation carbon emissions and time in a line graph.
EmissionsDropped.plot(subplots = True, figsize = (24,24), xlabel = 'Year', ylabel = 'Carbon Emissions (MMT)');
Observations:
- While coal may possibly have a constant mean and standard deviation, it has very strong seasonality which makes it non-stationary. The relatively constant mean indicates that coal has emitted roughly the same amount of carbon emissions over time.
- Distillate fuel may be stationary. There is no obvious seasonal pattern, the mean tends to stay constant when not spiking, and the standard deviation doesn't appear to change much. This may be wrong, depending on the impact of the sudden random spikes. Distillate fuel hasn't deviated much from about .25 MMT of carbon emissions every month over time.
- While geothermal's standard deviation appears constant and maybe the mean is constant, there is some strong seasonality. Overall, geothermal energy emits about .005 MMT more carbon emissions in 2016 than it did in 1989.
- Natural gas has some seasonality, a slightly increasing trend, and greater standard deviation with time. Natural gas emits about 53 MMT more carbon in 2016 than it did in 1989.
- Non-biomass waste has seasonality and a moving mean, but a constant standard deviation. The moving mean indicates an increase in carbon emissions of about .6 MMT.
- While petroleum coke doesn't appear to have a seasonal pattern, its standard deviation is greater with time and the mean appears to increase over time; it is emitting about .9 MMT more every month.
- Petroleum has no obvious seasonality, but the standard deviation and mean both decrease. It went from about 14 MMT a month in 1989 to about 1 MMT a month in 2016 for its carbon emissions.
- Residual fuel oil looks like petroleum no obvious seasonality and a decreasing standard deviation and mean. Residual fuel oil went from 12 MMT of CO2 in 1989 to about 1 MMT in 2016.
- Overall CO2 emissions appear to have a roughly constant mean around 175 MMT per month from 1989 to 2016. It does appear to have strong seasonality and the standard deviation increased slightly over time.
Rolling Averages and Standard Deviations Over Time¶
In [29]:
# Now we will view the rolling averages for all of the energy sources for the first 100 observations.
# Calculating the rolling mean and standard deviation for a window of 10 observations.
for c in EnergyTypes:
rolmean=EmissionsDropped[c].rolling(window=10).mean()
rolstd=EmissionsDropped[c].rolling(window=10).std()
#Visualizing the rolling mean and standard deviation
plt.figure(figsize=(16,8))
actual = plt.plot(EmissionsDropped[c], color='grey', label=c)
rollingmean = plt.plot(rolmean, color='blue', label='Rolling Mean')
rollingstd = plt.plot(rolstd, color='purple', label='Rolling Std. Dev.')
plt.title('Rolling Mean & Standard Deviation')
plt.legend()
plt.show()
Observations
- With seasonality in the rolling mean and a downward trend, coal is not a stationary series.
- Distillate fuel is difficult to tell if it is stationary due to the random spikes. The slight downward trend and changes in standard deviation indicate it is not a stationary series.
- Geothermal's strong seasonality keeps it from being a stationary series.
- Natural gas has an upward trend as well as a seasonal mean and standard deviation so it is not a stationary series.
- Non-biomass's upward trend and seasonality make it non-stationary.
- Petroleum Coke is probably not stationary due to the overall increase in its rolling mean.
- Petroleum's decreasing mean and slightly decreasing standard deviation keep it from being stationary.
- Residual oil is not stationary due to its decreasing standard deviation and rolling mean.
- Total emissions has seasonality in its rolling mean and a slight increase in standard deviation so it is not stationary.
Forcasting Preparation¶
For this project, we are only interested in analyzing natural gas CO2 emissions. Therefore, we need a dataframe which contains only the date as an index and the natural gas observations. When we went to visualize the data, we had to drop all of the observations measured before 1989 since some of the energy sources did not have data until then. However, natural gas has no missing values from the start of the original dataset, so we are going to use this to make our forcasting dataframe. Luckily, we still have observations from the entire daterange in our 'Emissions' dataframe!
In [30]:
# Create a new, natural gas-only dataframe based on the 'Emissions' dataframe.
NatGas = pd.DataFrame()
NatGas['Value'] = Emissions['Natural Gas']
NatGas.head(15)
Out[30]:
In [31]:
# View the dataframe's information.
NatGas.info()
Observations:
- This new dataframe has a single column of floats.
- There are 523 observations and no nulls in this series.
Milestone 1 Conclusion - Proposed Approach¶
- Potential techniques:
- To establish stationarity:
- Augmented Dickey-Fuller Test
- Logarithmic transformation
- Differencing the series
- Decompose the series into its trend, seasonality, and residual
- To model the shifted data:
- Use autocorrelation and partial correlation to determine p and q values
- Split the data 80/20 to create training and testing datasets, respectively
- Create the AR, MA, ARMA, and ARIMA models for the data
- Calculate the RMSE for each potential model
- Use inverse transformation to see how the model fares against the test data
- To establish stationarity:
- Overall solution design: The end solution should focus on only one of the models as determined by which one is the most successful.
- Measures of success:
- A low AIC value
- A low mean squared error
- No overfitting; a similar RMSE between the training and testing data.
Milestone 2¶
Splitting into training and testing data.¶
In [32]:
# Splitting the dataset into 80% train and 20% test and visualizing them on the same plot.
train, test = train_test_split(NatGas, test_size=.20, train_size=.80, random_state=10, shuffle=False);
fig, ax = plt.subplots(figsize = (24, 8));
train.plot(ax = ax, color = 'grey');
test.plot(ax = ax, color = 'purple');
plt.legend(['Train Data', 'Test Data']);
Establishing Stationarity¶
Next, we will implement the Augmented Dicky-Fuller Test to test for stationarity. The test has the following conditions:
- Null Hypothesis - the time series is non stationary due to seasonality and trend.
- Alternate Hypothesis - the time series is stationary.
- We are setting alpha at .05.
In [33]:
# We are going to visualize multiple plots of rolling means and standard deviations against the data.
# We will also use the Augmented Dicky-Fuller test more than once. To make this easier, we're going to make a function.
def MakePlot(x):
# Visualize the rolling mean and standard deviation with window = 12.
plt.figure(figsize = (24, 8))
plt.plot(x, color = 'gray', label = 'Values')
plt.plot(x.rolling(window = 12).mean(), color = 'blue', label = 'Mean', linestyle = '--')
plt.plot(x.rolling(window = 12).std(), color = 'purple', label = 'Std. Dev.', linestyle = ':')
plt.legend()
plt.show()
# Implement ADF test on the log data.
print(adfuller(x['Value']))
Testing for Stationarity¶
In [34]:
# Visualizing the rolling mean and standard deviation and performing the ADF test on the training data.
MakePlot(train)
Observations:
- The plot above shows strong seasonality with spikes mid-year and dips at the end and start of each year.
- It also shows that, overall, carbon emissions from natural gas have been gradually increasing since approximately 1990.
- The standard deviation is also increasing slightly over time.
- The series is not stationary so we need to decompose the time series into trend, seasonality, and white noise.
- Our p-value rounds up to 1 unless we include 6 decimal points. This is much higher than our alpha of .05 so we fail to reject the null hypothesis and say that the time series is not stationary.
In [35]:
train.info()
train.head()
Out[35]:
Visualise the decomposed data.¶
In [36]:
# Using seasonal_decompose to create a dataframe with the individual components of the decomposed time series.
decomposition = sm.tsa.seasonal_decompose(train.Value)
DecompData = pd.DataFrame()
DecompData['trend'] = decomposition.trend
DecompData['seasonal'] = decomposition.seasonal
DecompData['noise'] = decomposition.resid
# Visualize the components.
fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, figsize = (24,24));
DecompData['trend'].plot(ax=ax1, title='Trend', color='orange');
DecompData['seasonal'].plot(ax=ax2, title='Seasonal', color='purple');
DecompData['noise'].plot(ax=ax3, title='Noise');
Observations:
- There is an overall upward trend, indicating that there have been more natural gas carbon emissions over time.
- There is very strong seasonality with such a regular pattern of spiking in the summer and dropping in the winter that it looks intentional at a cursory glance.
- Noise appears to have a mean of 0, but its standard deviation may be increasing as time goes on.
- We need to use a combination of differencing and log transformation to create a stationary time series.
Perform logarithmic transformation¶
In [37]:
# We will use logarithmic transformation to stabilize the standard deviation and then test again for stationarity using ADF.
# Perform the logarithmic transformation.
LogTrain = np.log(train)
print(LogTrain)
# Visualizing the rolling mean and standard deviation and performing the ADF test on the transformed training data.
MakePlot(LogTrain)
Observations:
While the standard deviation has stabalized, there is still strong seasonality, an upward trend, and a p-value of .97. It is now time to shift the series by one month and apply differencing using the lagged series.
Perform first order differencing.¶
In [38]:
# Calculate the first order difference.
Diff = LogTrain.diff(periods = 1).dropna()
# Visualizing the rolling mean and standard deviation and performing the ADF test on the differenced training data.
MakePlot(Diff)
Observations:
The graph of the rolling mean and standard deviation are both stabalized and our p-value is only .001. Since this is less than our alpha of .05, we reject the null hypothesis and conclude the series is now stationary. We are ready to begin modeling the data.
In [39]:
# Using seasonal_decompose to create a dataframe with the individual components of the decomposed time series.
decomposition = sm.tsa.seasonal_decompose(Diff.Value)
DecompData = pd.DataFrame()
DecompData['trend'] = decomposition.trend
DecompData['seasonal'] = decomposition.seasonal
DecompData['noise'] = decomposition.resid
# Visualize the components.
fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, figsize = (48,24));
DecompData['trend'].plot(ax=ax1, title='Trend', color='orange');
DecompData['seasonal'].plot(ax=ax2, title='Seasonal', color='purple');
DecompData['noise'].plot(ax=ax3, title='Noise');
Observations:
There still appears to be some strong seasonality in the data, but our ADF p-value is already less than .05. Depending on how the AR, MA, ARMA, and ARIMA models perform, we may have to test other model types.
Build and Test Models¶
Find Optimal P and Q Parameters¶
In [40]:
# creating two subplots to show ACF and PACF plots
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20, 12))
# creating and plotting the ACF charts to lag=12
tsaplots.plot_acf(Diff, zero=False, ax=ax1, lags = 12)
# creating and plotting the ACF charts to lag=12
tsaplots.plot_pacf(Diff, zero=False, ax=ax2, lags = 12)
plt.show()
Observations:
- The ACF and PACF plots show that the highest lag at which each plot goes past the confidence interval boundary is lag 1.
- We therefore conclude that p = 1 and q = 1.
In [41]:
#Apply AR model
plt.figure(figsize=(16,8))
model_AR = AutoReg(Diff, lags=1) #Use number of lags as 1 and apply AutoReg function on Diff series
results_AR = model_AR.fit() #fit the model
plt.plot(Diff)
predict = results_AR.predict(start=0,end=len(Diff)-1) #predict the series
predict = predict.fillna(0) #Converting NaN values to 0
plt.plot(predict, color='red')
plt.title('AR Model - RMSE: %.4f'% mean_squared_error(predict,Diff['Value'], squared=False)) #Calculating rmse
plt.show()
In [42]:
results_AR.aic
Out[42]:
Observations:
So far we are off to a solid start with a low RMSE of 0.1539 and an AIC of about -369.51.
In [46]:
#Apply MA model.
plt.figure(figsize=(16,8))
model_MA = ARIMA(Diff, order=(0, 0, 1)) #Using p=0, d=0, q=1 and apply ARIMA function on Diff series
results_MA = model_MA.fit() #fit the model
plt.plot(Diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('MA Model - RMSE: %.4f'% mean_squared_error(results_MA.fittedvalues,Diff['Value'], squared=False))
plt.show()
In [47]:
results_MA.aic
Out[47]:
Observations:
While the MA model has both a higher RMSE and a higher AIC than the AR model.
In [48]:
# Apply the ARMA model.
plt.figure(figsize=(16,8))
model_ARMA = ARIMA(Diff, order=(1, 0, 1)) #Using p=1, d=0, q=1 and apply ARIMA function on Diff series
results_ARMA = model_ARMA.fit() #fit the model
plt.plot(Diff)
plt.plot(results_ARMA.fittedvalues, color='red')
plt.title('ARMA Model - RMSE: %.4f'% mean_squared_error(results_ARMA.fittedvalues,Diff['Value'], squared=False))
plt.show()
In [49]:
results_ARMA.aic
Out[49]:
Observations:
This model performed slightly worse than the AR model with an AIC of -369.32.
In [50]:
#Apply the ARIMA model.
plt.figure(figsize=(16,8))
model_ARIMA = ARIMA(Diff, order=(1, 1, 1))
results_ARIMA = model_ARIMA.fit() #fit the model
plt.plot(Diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('ARIMA Model - RMSE: %.4f'% mean_squared_error(results_ARIMA.fittedvalues,Diff['Value'], squared=False))
plt.show()
In [51]:
results_ARIMA.aic
Out[51]:
Observations:
ARIMA performed the worst of the for with a slightly higher RMSE and the highest AIC of -364.21. We will move forward with the AR model as it has the best RMSE and AIC.
In [66]:
# Printing the fitted values
predictions=pd.Series(results_AR.fittedvalues)
#First step - doing cumulative sum
predictions_cumsum = predictions.cumsum() # use .cumsum fuction on the predictions
#Second step - Adding the first value of the log series to the cumulative sum values
predictions_log = pd.Series()
predictions_log = pd.Series(LogTrain['Value'].iloc[13], index=LogTrain.index)
predictions_log = predictions_log.add(predictions_cumsum, fill_value=0)
#Third step - applying exponential transformation
predictions_AR = np.exp(predictions_log) #use exponential function
In [67]:
#Plotting the original vs predicted series
plt.figure(figsize=(16,8))
plt.plot(train, color = 'c', label = 'Original Series') #plot the original train series
plt.plot(predictions_AR, color = 'r', label = 'Predicted Series') #plot the predictions_ARIMA
plt.title('Actual vs Predicted')
plt.legend()
plt.show()
Observations:
While this model is doing a decent job of capturing the overall trend of the data, it is not accurately reflecting its seasonality. Let's see how it performs against the test data.
In [68]:
#Forecasting the values for the test data
forecasted_AR = pd.DataFrame()
forecasted_AR['forecasted'] = results_MA.forecast(steps=105)
#First step - doing cumulative sum
forecasted_cumsum = forecasted_AR.cumsum() # use .cumsum fuction on the predictions
#Second step - Adding the first value of the log series to the cumulative sum values
forecasted_log = pd.Series()
forecasted_log = pd.Series(LogTrain['Value'].iloc[417], index=forecasted_cumsum.index)
forecasted_log = forecasted_log + forecasted_cumsum['forecasted']
#Applying exponential transformation to the forecasted log values
forecasted_AR = np.exp(forecasted_log) #use exponential function on forecasted data
In [69]:
#Plotting the original vs predicted series
plt.figure(figsize=(16,8))
plt.plot(NatGas, color = 'c', label = 'Original Series')
plt.plot(predictions_AR, color = 'r', label = 'Prediction on Train data') #plot the predictions_ARIMA series
plt.plot(forecasted_AR, label = 'Prediction on Test data', color='b') #plot the forecasted_ARIMA series
plt.title('Actual vs Predicted')
plt.legend()
plt.show()
In [70]:
mean_squared_error(predictions_AR, train, squared = False) #calculate RMSE using the predictions_ARIMA and df_train
Out[70]:
In [71]:
mean_squared_error(forecasted_AR, test, squared = False) #calculate RMSE using the forecasted_ARIMA and df_test
Out[71]:
Observations:
The graph clearly demonstrates that our model failed to capture the seasonality of the data. This is confirmed by the fact that the forecasted RMSE is nearly double that of the training RMSE. Our model is overfitting.
Build and Test Sarima Model¶
Identify Optimal Parameters¶
We are going to brute force test for all combinations of all values of p, d, q, P, D, and Q between 0 and 2. While there may be better parameters out there, this has to be weighed against the time and resources cost of performing the test.
In [58]:
# Create a list of all possible parameter combinations.
p = d = q = range(0, 2)
Vals = list(itertools.product(p, d, q))
SarimaValues = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
# Perform logarithmic transformation on the original data.
Log = np.log(NatGas)
In [59]:
SarimaAics = []
for Param in Vals:
for SeasonalParam in SarimaValues:
try:
model_SARIMA = sm.tsa.statespace.SARIMAX(Log, order = Param, seasonal_order = SeasonalParam, enforce_stationarity=False, enforce_invertibility=False)
SarimaResults = model_SARIMA.fit()
print('ARIMA{}x{} - AIC:{}'.format(Param, SeasonalParam, SarimaResults.aic))
if SarimaResults.mle_retvals is not None and SarimaResults.mle_retvals['converged'] == False:
print(SarimaResults.mle_retvals)
SarimaAics.append(SarimaResults.aic)
except:
continue
SarimaAics.sort()
print(SarimaAics[0])
Observations:
Based purely on AIC, our best parameter combination is (1, 1, 1) x (1, 0, 1, 12) with an AIC of -1132.67. Before moving forward, though, we need to check and make sure that the residuals of the mean are uncorrelated and normally distributed with a zero-mean.
In [60]:
model_SARIMA = sm.tsa.statespace.SARIMAX(Log, order = (1, 1, 1), seasonal_order = (1, 0, 1, 12))
SarimaResults = model_SARIMA.fit()
print(SarimaResults.summary())
Observations:
- All of our model's terms have a statistically significant p-value, so we can keep all of them.
- The p-value of the Ljung-Box test is .50, which is not statistically significant. We therefore can't reject the null hypothesis that our residuals are white noise.
- The p-value of the heteroskedasticity is less than .05, so we reject the null hypothesis that the residuals are homoskedastic.
- The Jarque-Bera test's p-value is statistically significant, so we reject the null hypothesis that the residuals are normally distributed. This is due to the slight positive skew and large kurtosis.
In [61]:
SarimaResults.plot_diagnostics(figsize=(15, 12))
plt.show()
Observations:
- The standardized residual plot does not have any obvious seasonality and appears to be white noise.
- As the summary statistics indicated, the plot of the density function has a slight positive skew as the right tail has a local maximum.
- The Q-Q plot shows that the forecast errors slightly deviate from the red KDE line. This shows that the normal distribution is not a perfect fit for our residuals.
- All values in the correlogram are within the confidence interval.
- While the test statistics and plots show that our residuals do not perfectly resenble a normal distribution, it is not unreasonable. Given the excellent AIC of the SARIMA model and how slight the skew and kurtosis is for our residuals, the SARIMA model will probably give us better information and more accurate forecasts than the AR model would.
Testing the SARIMA Model.¶
Since we are going to test both the one-step-ahead forecast and the dynamic forecast as well as predict beyond the dates in our dataset, there are some processes we will perform repeatedly. To be more efficient, we will define functions to call later.
In [62]:
# This function gets predictions' confidence intervals, exponentiates them, plots them against the actual values, and prints the RMSE and MSE.
def PredResults(pred, errors = True):
# Get the fitted predictions.
predictions_SARIMAX = np.exp(pred.conf_int())
#Plot predictions agains values.
ax = NatGas['1973':].plot(label='observed',figsize=(24, 8))
np.exp(pred.predicted_mean).plot(label='Forecast', ax=ax, color = 'r')
ax.fill_between(predictions_SARIMAX.index,
predictions_SARIMAX.iloc[:, 0],
predictions_SARIMAX.iloc[:, 1], color='y', alpha=.5)
ax.set_xlabel('Time (years)')
ax.set_ylabel('NG CO2 Emissions')
plt.legend()
plt.show();
# Print the MSE and RMSE for predictions. If a forecast, check the p-values.
if errors == True:
NatGas_forecast = np.exp(pred.predicted_mean)
NatGas_truth = NatGas.Value.iloc[418:523]
mse = ((NatGas_forecast - NatGas_truth) ** 2).mean()
print('The Mean Squared Error (MSE) of the forecast is {}'.format(round(mse, 2)))
print('The Root Mean Square Error (RMSE) of the forecast: {:.4f}'
.format(np.sqrt(sum((NatGas_forecast-NatGas_truth)**2)/len(NatGas_forecast))))
else:
print(SarimaResults.pvalues)
return(pred, predictions_SARIMAX);
In [63]:
# Get predictions with a one-step-ahead forecast.
OneStep = PredResults(SarimaResults.get_prediction(start = 418, end = 522, dynamic=False))
Observations:
While the RMSE is higher for the SARIMA model than it was for the AR model, we can clearly see from the graph that the SARIMA model did a much better job capturing the magnitude of the data's seasonality. Maybe a dynamic model will perform better.
In [64]:
# Get predictions with a dynamic forecast.
Dynamic = PredResults(SarimaResults.get_prediction(start = 418, end = 522, dynamic=True))
Observations:
The dynamic SARIMA's MSE and RMSE are 4 times and 2 times greater than those of the one-step-ahead forecast, respectively. We will therefore use the one-step-ahead forecast to predict the next 2 years.
In [65]:
# Forecast the next two years with a one-step-ahead forecast.
Forecast = PredResults(SarimaResults.get_forecast(steps = 24, dynamic=False), False)
Module 2 Conclusion¶
Final Solution Design:
- We have tested the CO2 emissions forecasting accuracy of the AR, MA, ARMA, ARIMA, and SARIMA methods.
- SARIMA's low RMSE, very low AIC, fulfillment of all necessary statistical assumptions, and ability to model the data's seasonality make it the best model out of those tested.
- SARIMA is the preferred model. The RMSE and the Jaque-Bera test both suggest that more refinement beyond the scope of this project could result in an even better model.
- We have forecasted the next 24 months of CO2 emissions data.
Refined Insights:
- Natural gas CO2 emissions have a very strong annual seasonality.
- Coal currently emits more CO2 than any other energy source, including natural gas.
- This may soon change as coal's CO2 emissions started declining around 2008 and natural gas's CO2 emissions began increasing around 2000. If coal were to continue decreasing from its CO2 emissions of approximately 100 MMT and the forecast for natural gas is accurate, it is highly possible natural gas will become the greatest CO2 contributor.
- We failed to reject the null hypothesis that the series is stationary due to its variation in standard deviation and mean over time. We therefore need to be cautious and not assume that natural gas's emissions are dependent on coal's when it could merely be dependent on time.
- It would behoove future researchers to take into account the average Watt to ton ratio of each of the energy sources and who is using which. It could be better or worse for the environment that natural gas appears to be slowly replacing coal.